import os
import numpy as np
import pandas as pd
import scipy.stats as sps
from tqdm.auto import tqdm
from sklearn.decomposition import PCA

import matplotlib.pyplot as plt
import seaborn as sns
import cv2
import imageio
from IPython.display import Image

sns.set_theme(style="darkgrid", palette="Set2")

# Сгенерируем 100 точек в двумерном пространстве
n = 100  # Размер выборки
sample = sps.multivariate_normal(mean=[5, -7], cov=[[13, 7], [7, 5]]).rvs(
    n, random_state=42
)

# Визуализируем полученную выборку
plt.figure(figsize=(10, 6))
sns.scatterplot(x=sample[:, 0], y=sample[:, 1], s=50, c="g")

plt.title("Двумерные синтетические данные")
plt.xlabel("Первая размерность")
plt.ylabel("Вторая размерность")
plt.show()

pca = PCA(n_components=1)
# Обучаем и тут же применяем PCA к выборке sample
sample_projection = pca.fit_transform(sample)

print(
    f"Направления главных компонент в исходном пространстве (одной компоненты): {pca.components_.round(3)}"
)
print(f"Среднее по данным: {pca.mean_.round(3)} == {sample.mean(axis=0).round(3)}")
print(f"Размерность проекции: {sample_projection.shape}")

# Визуализируем на числовой прямой проекцию после применения PCA
plt.figure(figsize=(10, 1))
sns.scatterplot(x=sample_projection[:, 0], y=0, c="b", alpha=0.5)

plt.title("Проекция исходных данных на первую главную компоненту")
plt.show()

Направления главных компонент в исходном пространстве (одной компоненты): [[0.87  0.494]]
Среднее по данным: [ 5.396 -6.732] == [ 5.396 -6.732]
Размерность проекции: (100, 1)

# Применяем обратное преобразование
sample_reprojection = pca.inverse_transform(sample_projection)
print(f"Размерность репроекции в исходное пространство: {sample_reprojection.shape}")

# Визуализируем
plt.figure(figsize=(10, 6))

# Рисуем отрезки, соединяющие исходную точку и ее образ
for point, reproj in zip(sample, sample_reprojection):
    plt.plot([point[0], reproj[0]], [point[1], reproj[1]], c="r", alpha=0.3)

# Рисуем исходные точки
sns.scatterplot(x=sample[:, 0], y=sample[:, 1], s=50, c="g", label="Исходные данные")
# Рисуем восстановленные точки
sns.scatterplot(
    x=sample_reprojection[:, 0],
    y=sample_reprojection[:, 1],
    s=25,
    c="b",
    label="Репроекция после PCA",
)

plt.axis(
    "equal"
)  # Включаем одинаковый масштаб по обоим осям, чтобы увидеть прямой угол
plt.title("Исходные данные и репроекция после PCA")
plt.xlabel("Первая размерность")
plt.ylabel("Вторая размерность")
plt.show()

Размерность репроекции в исходное пространство: (100, 2)

faces = np.loadtxt("data/train_faces.npy")
print(f"Размерность faces: {faces.shape}")

# Визуализируем несколько примеров из датасета
plt.figure(figsize=(12, 5))
for i in range(40):
    plt.subplot(4, 10, i + 1)
    plt.imshow(faces[i].reshape((32, 32)), cmap="gray")
    plt.axis("off")

Размерность faces: (280, 1024)

pca = PCA(n_components=50)
faces_projection = pca.fit_transform(faces)
print(f"Размерность проекции: {faces_projection.shape}")

Размерность проекции: (280, 50)

def rescale_image(image):
    """
    Линейно переводит минимальное значение массива в 0, а максимальное в 1

    :param image: массив для преобразования
    """
    return (image - image.min()) / (image.max() - image.min())


def visualize_main_components(pca, num_components=50, image_shape=(32, 32)):
    """
    Визуализирует первые главные компоненты при PCA проекции

    :param pca: обученный PCA
    :param num_components: сколько компонент визуализировать
    :param image_shape: размерность изображения, например,
                        (32, 32) для ч/б или (64, 64, 3) для цветного изображения
    """
    cmap = "gray" if image_shape[-1] != 3 else None
    n_cols = 10
    n_rows = (num_components + n_cols - 1) // n_cols

    plt.figure(figsize=(1, 1))
    plt.imshow(rescale_image(pca.mean_.reshape(image_shape)), cmap=cmap)
    plt.title("Среднее")
    plt.axis("off")
    plt.show()

    plt.figure(figsize=(n_cols * 1.1, n_rows * 1.1))
    plt.suptitle(
        f"Первые {num_components} главных компонент в порядке убывания важности"
    )
    for i in range(num_components):
        plt.subplot(n_rows, n_cols, i + 1)
        # Так как значения каждого пикселя изображения должны быть на отрезке [0, 1],
        # а главные компоненты могут не удовлетворять этому условию, отнормируем их
        plt.imshow(rescale_image(pca.components_[i].reshape(image_shape)), cmap=cmap)
        plt.axis("off")

    plt.tight_layout()


visualize_main_components(pca, num_components=50)

def visualize_projection_images(
    images, n_components, num_images=5, image_shape=(32, 32)
):
    """
    Визуализирует проекции изображений на первые несколько главных компонент

    :param images: датасет изображений, на котором обучается PCA
    :param n_components: количество главных компонент
    :param num_images: количество первых изображения датасета для отображения
    :param image_shape: размерность изображения
    """

    pca = PCA(n_components=n_components)
    projection = pca.fit_transform(images)
    reprojection = pca.inverse_transform(projection)

    cmap = "gray" if image_shape[-1] != 3 else None
    num_columns = 3 + n_components

    plt.figure(figsize=(1.1 * num_columns, 1.2 * num_images))
    plt.suptitle(
        f"Веса при разложении изображений на первые {n_components} главных компонент"
    )

    for image_index in range(num_images):
        plt.subplot(num_images, num_columns, image_index * num_columns + 1)
        plt.imshow(images[image_index].reshape(image_shape), cmap=cmap)
        plt.title("лицо")
        plt.axis("off")

        plt.subplot(num_images, num_columns, image_index * num_columns + 2)
        plt.imshow(
            np.clip(reprojection[image_index].reshape(image_shape) / 255, 0, 1),
            cmap=cmap,
        )
        plt.title("прибл.")
        plt.axis("off")

        plt.subplot(num_images, num_columns, image_index * num_columns + 3)
        plt.imshow(np.clip(pca.mean_.reshape(image_shape) / 255, 0, 1), cmap=cmap)
        plt.title("= 1")
        plt.axis("off")

        for component_idx in range(n_components):
            plt.subplot(
                num_images, num_columns, image_index * num_columns + 4 + component_idx
            )
            plt.imshow(
                rescale_image(pca.components_[component_idx].reshape(image_shape)),
                cmap=cmap,
            )

            coef = projection[image_index, component_idx]
            sign = "+" if coef >= 0 else ""
            plt.title(f"{sign}{coef:.1f}")
            plt.axis("off")

    plt.tight_layout()
    plt.show()


visualize_projection_images(faces, n_components=10, num_images=10, image_shape=(32, 32))

def visualize_decompress_images(
    images, n_components_grid, num_images=5, image_shape=(32, 32)
):
    """
    Визуализирует восстановленные изображения для различного количества главных компонент

    :param images: датасет изображений, на котором обучается PCA
    :param n_components_grid: сетка количества главных компонент
    :param num_images: количество первых изображения датасета для отображения
    :param image_shape: размерность изображения
    """
    reprojections = {}

    for n_components in tqdm(n_components_grid):
        pca = PCA(n_components=n_components)
        projection = pca.fit_transform(images)
        reprojections[n_components] = pca.inverse_transform(projection)

    num_columns = 2 + len(n_components_grid)

    cmap = "gray" if image_shape[-1] != 3 else None
    plt.figure(figsize=(1.2 * len(n_components_grid), 1.1 * num_images))
    plt.suptitle(
        "Зависимость восстановленных изображений от количества главных компонент"
    )

    for image_index in range(num_images):
        plt.subplot(num_images, num_columns, image_index * num_columns + 1)
        plt.imshow(images[image_index].reshape(image_shape), cmap=cmap)
        if image_index == 0:
            plt.title("лицо")
        plt.axis("off")

        plt.subplot(num_images, num_columns, image_index * num_columns + 2)
        plt.imshow(np.clip(pca.mean_.reshape(image_shape) / 255, 0, 1), cmap=cmap)
        if image_index == 0:
            plt.title("среднее")
        plt.axis("off")

        for col, n_components in enumerate(n_components_grid):
            plt.subplot(num_images, num_columns, image_index * num_columns + 3 + col)
            reproj = reprojections[n_components][image_index]
            plt.imshow(np.clip(reproj.reshape(image_shape) / 255, 0, 1), cmap=cmap)
            if image_index == 0:
                plt.title(f"{n_components} комп.")
            plt.axis("off")

    plt.tight_layout()
    plt.show()

    return reprojections


faces_reprojections = visualize_decompress_images(
    faces,
    n_components_grid=[1, 5, 10, 20, 50, 75, 100, 150, 200],
    num_images=10,
    image_shape=(32, 32),
)

cats_path = "data/cats"
cats = []

for filename in os.listdir(cats_path):
    cats.append(plt.imread(os.path.join(cats_path, filename)))

cats = np.stack(cats).reshape(-1, 64 * 64 * 3)
print(f"Размерность cats: {cats.shape}")

# Визуализируем примеры
plt.figure(figsize=(12, 5))
for i in range(40):
    plt.subplot(4, 10, i + 1)
    plt.imshow(cats[i].reshape(64, 64, 3))
    plt.axis("off")

Размерность cats: (15747, 12288)

pca = PCA(n_components=100)
cats_projection = pca.fit_transform(cats.astype(np.float32))
print(f"Размерность проекции: {cats_projection.shape}")

visualize_main_components(pca, num_components=70, image_shape=(64, 64, 3))

Размерность проекции: (15747, 100)

visualize_projection_images(
    cats, n_components=10, num_images=10, image_shape=(64, 64, 3)
)

cats_reprojections = visualize_decompress_images(
    cats,
    n_components_grid=[1, 5, 10, 20, 50, 100, 200, 500, 1000],
    num_images=10,
    image_shape=(64, 64, 3),
)

pca = PCA(n_components=1000)
cats_projection = pca.fit_transform(cats.astype(np.float32))

cnt_cats = 9
gif_images = []

image_size = 64
border_size = 2

for num_components in range(1, 1000):
    # Делаем проекцию на первые num_components главных компонент
    projection = (cats[:cnt_cats] - pca.mean_) @ pca.components_[:num_components].T
    # Восстанавливаем спроецированные изображения
    reprojection = projection @ pca.components_[:num_components] + pca.mean_

    # Решейпим вектора в матрицы пикселей
    compressed_images = (
        reprojection.reshape(-1, image_size, image_size, 3)
        .clip(0, 255)
        .astype(np.uint8)
    )

    # Подготавливаем главную компоненту
    pca_component = (
        rescale_image(pca.components_[num_components].reshape(64, 64, 3)) * 255
    ).astype(np.uint8)

    # Собираем коллаж 3х3: в центре текущая num_components-ая главная компонента,
    # по бокам восстановленные изображения
    collage = np.zeros(
        (image_size * 3 + border_size * 2, image_size * 3 + border_size * 2, 3),
        dtype=np.uint8,
    )
    for i, image in enumerate(compressed_images):
        row_begin = (i % 3) * (border_size + image_size)
        col_begin = (i // 3) * (border_size + image_size)

        # Если i == 4, то это центральное изображение, заменяем его на главную компоненту
        if i == 4:
            image = pca_component
        collage[
            row_begin : row_begin + image_size, col_begin : col_begin + image_size
        ] = image

    # Добавляем полученный коллаж в общий массив, перед этим растянув его в два раза
    gif_images.append(
        cv2.resize(
            collage,
            dsize=(collage.shape[0] * 2, collage.shape[0] * 2),
            interpolation=cv2.INTER_CUBIC,
        )
    )

# Для ускорения гифки на большом количестве компонент будем выкидывать часть кадров
frames = gif_images[:50] + gif_images[50:200:3] + gif_images[200::10]
imageio.mimsave("transformation.gif", frames * 5, format="GIF")

Image(url="transformation.gif")

Введение в анализ данных ¶

Метод главных компонент (Principal Component Analysis, PCA)¶

Реализация PCA из sklearn¶

1. Простой пример¶

2. PCA для сжатия датасетов изображений лиц¶

3. PCA для сжатия датасетов изображений котиков¶

4. Бонус: гифка¶

Введение в анализ данных¶

Метод главных компонент (Principal Component Analysis, PCA)¶

Реализация PCA из sklearn¶

1. Простой пример¶

2. PCA для сжатия датасетов изображений лиц¶

3. PCA для сжатия датасетов изображений котиков¶

4. Бонус: гифка¶

Введение в анализ данных ¶